home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_9 / a9_6.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  3.2 KB  |  138 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 9.6 (Adams-Bashforth-Moulton Method).
  9. % Section    9.6,    Predictor-Corrector Method, Page 471
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program implements Adams-Bashforth-Moulton method
  13.  
  14. % for solving the initial value problem
  15.  
  16. %     y' = f(t,y)    with    y(a) = y0
  17.  
  18. %    Define and store the function f(t,y)    in the M-file  f.m 
  19. % function z = f(t,y)
  20. % z = 30 - 5*y;
  21.  
  22. delete f.m
  23. diary f.m; disp('function z = f(t,y)');...
  24.            disp('z = 30 - 5*y;');...
  25. diary off;
  26.  
  27. f(0,0); % Remark.  f.m rk4.m abm.m are used for Algorithm 9.6
  28. pause      % Press any key to continue.
  29.  
  30. clc;
  31. %    Place the endpoints of [a,b] in  a  and  b.
  32.  
  33. %    Place the initial value ya = y(a) in ya.
  34.  
  35. %    Place the number of subintervals in  n.
  36.  
  37. a  = 0;
  38. b  = 10;
  39. ya = 1;
  40. n  = 65;
  41.  
  42. pause    % Press any key to continue. 
  43.  
  44. clc;
  45. % Set up the step size, and vectors T and Y.
  46. h  = (b - a)/n;
  47. T = zeros(1,n+1);
  48. Y = zeros(1,n+1);
  49.  
  50. % Generate three starting values using the Runge-Kutta method.
  51. [Ts,Ys] = rk4('f',a,a+3*h,ya,6);
  52. T(1:4) = Ts(1:2:7);
  53. Y(1:4) = Ys(1:2:7);
  54.  
  55. % Proceed with the Adams-Bashforth-Moulton method.
  56. [T,Y]  = abm('f',T,Y);
  57. points = [T;Y];
  58.  
  59. pause    % Press any key to see the list of solution points.
  60.  
  61. clc; clg;
  62. c = min(Y); c = 0; % Set range.
  63. d = max(Y); d = 7; % Set range.
  64. axis([a b c d]);...
  65. plot(T,Y,'-r');...
  66. hold on;...
  67. plot([a b],[0 0],'b',[0 0],[c d],'b');...
  68. xlabel('t');...
  69. ylabel('y');...
  70. Mx1 = 'Adams-Bashforth-Moulton solution to y` = f(t,y).';...
  71. title(Mx1);...
  72. grid;...
  73. shg; pause    % Press any key to continue.
  74.  
  75. Mx2 = '     t(k)               y(k)';
  76. Mx3 = 'The Adams-Bashforth-Moulton method is stable';
  77. Mx4 = ['for  n = ',num2str(n),'  and  h = '];
  78. Mx5 = [Mx4,num2str(h),' because'];
  79. Mx6 = '     h < 0.75/|f (t,y)|';
  80. Mx7 = '                y';
  81. clc,echo off,diary output,...
  82. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  83. disp(Mx3),disp(Mx5),disp(''),disp(Mx6),disp(Mx7),...
  84. diary off, echo on
  85.  
  86. pause    % Press any key to continue.
  87.  
  88. clc;
  89. %    Place the endpoints of [a,b] in  a  and  b.
  90.  
  91. %    Place the initial value y(a) in ya.
  92.  
  93. %    Place the number of subintervals in  n.
  94.  
  95. a  = 0;
  96. b  = 10;
  97. ya = 1;
  98. n  = 37;
  99.  
  100. pause    % Press any key to continue. 
  101.  
  102. clc;
  103. % Set up the step size, and vectors T and Y.
  104. h  = (b - a)/n;
  105. T = zeros(1,n+1);
  106. Y = zeros(1,n+1);
  107.  
  108. % Generate three starting values using the Runge-Kutta method.
  109. [Ts,Ys] = rk4('f',a,a+3*h,ya,6);
  110. T(1:4) = Ts(1:2:7);
  111. Y(1:4) = Ys(1:2:7);
  112.  
  113. % Proceed with the Adams-Bashforth-Moulton method.
  114. [T,Y]  = abm('f',T,Y);
  115. points = [T;Y];
  116.  
  117. pause    % Press any key to see the list of solution points.
  118.  
  119. clc;
  120. plot(T,Y,'-g');...
  121. xlabel('t');...
  122. ylabel('y');...
  123. title(Mx1);...
  124. grid;...
  125. axis;...
  126. hold off;...
  127. shg; pause    % Press any key to continue.
  128.  
  129. Mx3 = 'The Adams-Bashforth-Moulton method is unstable';
  130. Mx4 = ['for  n = ',num2str(n),'  and  h = '];
  131. Mx5 = [Mx4,num2str(h),' because'];
  132. Mx6 = '     0.75/|f (t,y)| < h';
  133. Mx7 = '            y';
  134. clc,echo off,diary output,...
  135. disp(''),disp(Mx1),disp(''),disp(Mx2),disp(points'),...
  136. disp(Mx3),disp(Mx5),disp(''),disp(Mx6),disp(Mx7),...
  137. diary off, echo on
  138.